## 

rm(list = ls())

## 

library(tidyverse)
library(cobalt)
library(knitr)
library(kableExtra)
library(rdrobust)

## Get data

bt <- readRDS("data/data_federal.rds") %>%
  filter(!is.na(treated)) %>%
  mutate(year = lubridate::year(date.id)) %>%
  mutate(state_id = substr(ags, 1, 2))

## Drop 2009

bt <- bt %>% 
  filter(year > 2012)

## Covariate list

cov_list <- c('pop_foreign_share',
              'pop_share_65plus_2011',
              'soz_vers_beschaeftigte_share',
              'unem_capita',
              'households_prop_married',
              'households_size_2ormore',
              'residential_share_owneroccupied_2011', 
              'residential_share_100sqmplus_2011', 
              'residential_share_2000later_2011', 
              'pop_density_km2', 
              'migration_out_share', 
              'right_total_party_2009_btw',
              'left_total_party_2009_btw',
              'turnout_party_2009_btw', 
              'euro_per_sqm_county_2012', 
              'gdp_capita_county_2012')

## Source helper function to tidy rdrobust output

source("code/tidy_rd.R")

## Prep for RD

bal_df <- bt %>% 
  filter(!state_id == '08') %>% 
  filter(year == 2017) %>%
  dplyr::select(ags, pop_dec_09, applies_census,
                one_of(cov_list), year) %>%
  distinct(ags, .keep_all = T) %>%
  mutate(runvar = (pop_dec_09 * -1) + 10000) %>%
  filter(applies_census == 1) %>% ungroup()

## Standardize everything

bal_df <- bal_df %>% 
  mutate_at(vars(one_of(cov_list)), list(~(. / sd(., na.rm = T))))

## Add treatment

bal_df <- bal_df %>% 
  mutate(treated = ifelse(pop_dec_09 < 10000, 1, 0))

## Get RD estimates for the covariates

bw_ols = 2500

out_rd_opt <- pblapply(cov_list, function(o) {
  
  ## Estimate via RD
  
  out <- rdrobust(y = bal_df[, o] %>% pull(!!o), 
                  x = bal_df$runvar, c = 0)
  
  ## Tidy and return
  
  out2 <- out %>% 
    tidy_rd(se_nr = 3) %>% 
    mutate(outcome = o, method = 'RDD (MSE-optimal BW)') %>% 
    dplyr::rename(se = std.error, pval = p.value, 
                  bw = bw_left_h,
                  bw_bias = bw_left_b) %>% 
    mutate(conf.low90 = estimate - qnorm(0.95)*se,
           conf.high90 = estimate + qnorm(0.95)*se)
  
  ## Both
  
  out2
}) %>%
  reduce(rbind) 

## Prep data for table

## Proper labels

new_cov_list <- c('Foreign born / capita (2011)',
                  'Age 65+ / capita (2011)',
                  'Employment / capita (2011)',
                  'Unemployment / capita (2011)', 
                  'Households: prop. married couples (2011)',
                  'Households: prop. 2+ members (2011)',
                  'Residences: prop. owner-occupied (2011)',
                  'Residences: prop. 100+ sqm area (2011)',
                  'Residences: prop. built 2000 or later (2011)',
                  'Population density / km2 (2011)',
                  'Out-migration / capita (2011)', 
                  'Right party vote share (2009)', 
                  'Left party vote share (2009)', 
                  'Turnout (2009)',
                  'Land value (county-level, 2012)',
                  'GDP / capita (county-level, 2012)')

## Rename

for (i in 1:length(cov_list)) {
  out_rd_opt$outcome[out_rd_opt$outcome == cov_list[i]] <- 
    new_cov_list[i]
}


## Make table df

table_df <- out_rd_opt %>% 
  mutate_at(vars(estimate, matches('conf'), 
                 pval), list(~round(., 3))) %>% 
  mutate(ci = paste0('[',conf.low, ', ',conf.high, ']')) %>% 
  dplyr::select(outcome, method, estimate, 
                ci, pval, bw, bw_bias, n) %>%
  mutate(bw = round(bw, 0), bw_bias = round(bw_bias, 0))

## Table 2: Balance ----

kable(table_df %>% 
        filter(!str_detect(method, 'OLS')) %>% 
        dplyr::select(-method), "latex", longtable = F, 
      booktabs = T, col.names = c('Outcome',
                                  '$\\hat{\\tau}_{\\text{SRD}}$',
                                  'CI', 
                                  'P-val', 
                                  '$h_{\\text{MSE}}$',
                                  '$b_{\\text{MSE}}$',
                                  '$n$'),
      linesep = "",
      escape = F) %>%
  kable_styling(latex_options = c("repeat_header"),
                font_size = 9) %>% 
  # collapse_rows() %>% 
  row_spec(0, bold = T) %>% 
  kable_styling(latex_options = "HOLD_position")

